BI-TP 2012/13 



Charmonium properties in hot quenched lattice QCD 

H.-T. A. Francis^-^, O. Kaczmarek^ F. Karsch^'^^ 

H. Satzi and W. Socldner^ 
^Fakultdt fiir Physik, Universitdt Bielefeld, D-33615 Bielefeld, Germany 
^Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA 
^Institut fiir Kernphysik, Johannes Gutenberg Universitdt Mainz, D- 55099 Mainz, Germany 
"^Institut fiir Theoretische Physik, Universitdt Regensburg, D-93040 Regensburg, Germany 

We study the properties of charmonium states at finite temperature in quenched QCD on large 
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functions both below and above Tc. Our analysis suggests that both S wave states (J/V' and 
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and approximately I/ttT at 1.5 < T < 3 T^. 
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I. INTRODUCTION 



The main goal of ongoing heavy ion programs at the RHIC and LHC is to study the properties of the hot and dense 
medium formed during the colhsion of two relativistic heavy nuclei. It has been conjectured that at sufficiently high 
temperature the QCD medium will undergo a phase transition to a deconfined phase, in which the degrees of freedom 
are those of quarks and gluons. Unlike light mesons, the heavy mesons, e.g. '/A/', may survive in the hot medium 
up to certain temperatures before they get dissociated due to Debye screening Thus due to the dissociation of 
heavy mesons, the suppression of their yield in nucleus-nucleus (AA) collisions compared to that in proton-proton (pp) 
collisions can serve as a good probe for the properties of the medium. The experiments carried out at the SPS and 
LHC at CERN and the RHIC at BNL have indeed observed J/ijj suppression [31 . The interpretation of experimental 
data, however, is not as straightforward as the original idea proposed in Ref. The observed modification of J/ip 
production in AA collisions could be caused by two distinct classes of effects. On the one hand there are cold nuclear 
matter effects, which originate from the presence of cold nuclear matter in the target and projectile. On the other 
hand there are hot medium effects, which are of primary interest and reflect the properties of the medium we want 
to study. In order to disentangle these two effects, it is crucial to have a good understanding of the behavior of heavy 
quarks and quarkonia in the hot medium. 

From the theoretical point of view, the meson spectral function at finite temperature Q, which contains all the 
information on the hadron properties in the thermal medium, such as the presence, the location and the width of 
bound states (and thus about dissociation temperatures) as well as transport properties (e.g. heavy quark diffusion 
coefficients) , is the key quantity to be investigated. As this is a difficult task, several theoretical approaches have been 
followed to determine the quarkonium properties at finite temperature. 

The most traditional approach is based on the analysis of nonrelativistic potential models. Here one assumes that 
the interaction between a heavy quark pair forming the quarkonium can be described by a potential Q. Because 
of its success at zero temperature, the potential model approach has been used also at finite temperature @. The 
temperature dependent potential used in these calculations is based either on model calculations or on finite temper- 
ature lattice QCD results It is used to solve a nonrelativistic Schrodinger equation. The resulting dissociation 
temperatures depend strongly on the potential used. Recently progress has been made in comparing directly heavy 
quark correlation functions calculated on the lattice with potential model results. This allows to eliminate certain 
ambiguities [Bj and opens the possibility to determine which potential is more appropriate for a description of the 
experimental data 0. Nonetheless, the potential model approach at finite temperature is still under scrutiny. 

Most recently a nonrelativistic effective theory approach at nonzero temperature, which requires the scales concerned 
to be in hierarchy, has been developed Q. By integrating out certain scales, one arrives at a complex real-time static 
potential, which includes effects of screening via its real part as well as the interaction with the medium via its 
imaginary part. The presence of an imaginary part in the heavy quark potential reduces the possibility for stable 
quarkonium states in the hot medium. This approach becomes more reliable as the quark mass increases and thus is 
more relevant for the analysis of bottomonium states. An approach to study charmonium spectral functions at finite 
temperature using QCD sum rules has also been developed recently Q. 

First principle calculations in lattice QCD are thus crucially needed to determine the nonperturbative behavior of 
heavy quarks and quarkonia in the hot medium. The investigations of charmonium states at finite temperature, which 
have been performed in both quenched and full lattice QCD, have led to the rather interesting result that J/ip appears 
to survive up to temperatures well above Tc [l0l - [l6| . The most relevant quantities, meson spectral functions, however, 
cannot be obtained directly from lattice QCD calculations. Further input is needed to extract spectral functions 
from correlation functions calculated on the lattice. One of the commonly used methods is the Maximum Entropy 
Method (MEM). When using the Maximum Entropy Method, a very important issue is to get control over its input 
parameter (default model) dependence. The output spectral function from MEM can only be trusted if the default 
model dependence is eliminated or at least well understood. Additionally sufficient information on the Euclidean time 
dependence of the correlation function is crucially important in the MEM analysis. One economical way to increase 
the number of correlator data points in the temporal direction is to perform simulations on anisotropic lattices [TTI - 
[l6j . However, lattice cutoff effects are more significant on anisotropic lattices (l7| . Thus in the present paper we use 
isotropic lattices and perform simulations on very large lattices. The finest lattices we performed simulations on are 
128^ X 96, 128^ X 48, 128^ x 32 and 128^ x 24 at 0.73 T^, 1.46 T^, 2.20 and 2.93 T^, respectively. The number 
of data points in the temporal direction in the current study is doubled compared to our previous study in Ref. [lo| 



and it is about 1.5 times larger than that in stuidies [lll - ll5| and compatible with that used in Ref.[16|. Based on 
the correlation functions calculated on these large lattices, we will report on a detailed study of finite temperature 
charmonium correlators and perform a detailed MEM analysis of spectral functions, expanding on preliminary results 
reported in [l8l - [2l| . The signature obtained for the dissociation of charmonium states in the hot medium from the 
spectral function will be discussed. 

Besides the properties of charmonium states in the medium, the behavior of a single charm quark in the hot medium 





3 



is of great interest as well. Experimentally a substantial elliptic flow of heavy quarks has been observed The 
heavy quark diffusion D can be connected to the energy loss of a heavy quark during its propagation in the medium 
and is also related to the ratio of shear viscosity to entropy density ry/s [1^ Various phenomenological model 
studies suggest the heavy quark diffusion coefficient D < 1/T to accommodate data while various pQCD and T-Matrix 
calculations of the heavy quark diffusion coefficient differ significantly from each other [23 - [27| . It is thus important 
to have a first principle calculation of the heavy quark diffusion coefficient. The heavy quark diffusion coefficient can 
be obtained from the vector spectral function at vanishing frequency through the Kubo formula. We will give here 
also an estimate for the value of the charm diffusion coefficient at different temperatures. 

The rest of the paper is organized as follows. In Sec. |lT]we discuss general features of quarkonium correlators and 
spectral functions. In Sec. lIIII we give the lattice setup used in the calculation of charmonium correlation functions. In 
Sec. II VI we discuss information on the change of spectral functions from below to above that can be obtained from 
the analysis of correlation functions only, i.e. on thermal modifications of charmonium states and also on the charm 
quark diffusion coefficient. In Sec. |V]we will describe the Maximum Entropy Method used for the reconstruction 
of spectral functions and discuss the spectral functions below and above Tc obtained from MEM. Signatures for the 
dissociation of charmonium states and values of charm quark diffusion coefficients are discussed. Finally we summarize 
in Sec. IVII Some further details of our MEM analyses are given in an Appendix. 



II. MESON CORRELATION AND SPECTRAL FUNCTIONS 



In this section, we give the definition of the meson spectral function and its relation to the Euclidean correlation 
function, which can be calculated directly on the lattice. 

All information on quarkonium states is embedded in these spectral functions. The spectral function for a given 
meson channel H in a system can be defined through the Fourier transform of the real-time two-point correlation 
functions and D~ . The ensemble average of the commutator is 

Dnit, x) = -i {[Jnit, x), J^(0, 0)]) ^ D+{t, x) - D]j{t, x), (I) 

and its spectral density p{uj, p) can be expressed in terms of the retarded correlator D^{uj, p) [2^ 

PH{io,p) ^ D+{io,p) - D],{io,p) =2lmD^{io,p), (2) 

where 

Dp~\co,p)^ Jd^xe^"'-^P-Dp-\t,^), (3) 

i£)+(t,x) = (Jff(t,x)JH(0,O)), (4) 
ii^^(i,x) = (JH(0,O)JH(t,x)) . (5) 

The two-point correlation functions D^'" satisfy the Kubo-Martin-Schwinger (KMS) relation 

i^+(t,x) = i?^(t + */3,x), D+iu:,p)^e^^D],{Lu,p). (6) 
Inserting a complete set of states in Eq. ^ and using the KMS relation, one gets an explicit expression for ph{^, p) 

p„{cj,p) = ^^e"^^" (^6{p + kn-k,n)-6{p + k,n-kn)) | (n | (0) | m) | ' , (7) 



n.7n 



where Z is the partition function, p = (a;,p) and fcn(m) refer to the four-momenta of the state |n(m)). Given the 
above equation it is clear that the spectral function ph{uj,p) is an odd function of the frequency and momentum, 
Ph{—uj, — p) = —pfj{ijj, p) and ujpfj{uj, p) > 0. If the system is rotationally invariant, which means the state can have 
the same energy a; but opposite momentum p, the spectral function /0//(w, p) would also be an odd function of lu. 

The spectral function in the vector channel is related to the experimentally accessible differential cross section for 
thermal dilepton production (29| . 

dW 5a2 1 

■pv(w,p,r), (8) 



dcjdSp 547r3a;2(e"/^ - 1) 

where a is the electromagnetic fine structure constant and pv is the spectral function in the vector channel. Addi- 
tionally the spatial components of the vector spectral function are related to the heavy quark diffusion constant D Q 
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i=i 

where x^'^ is the quark number susceptibihty that is defined through the zeroth component of the temporal correlator 
in the vector channel. 

In this work we consider local meson operators of the form 

Jh(t,x) =^(T,x)rH^(T,x), (10) 

with Th = 1,75,7m,757ai, for scalar {SC), pseudoscalar {PS), vector {VC) and axial- vector {AV) channels, re- 
spectively. The relation of these quantum numbers to different charmonium states from the particle data book is 
summarized in Table HI 



Channel 


Th 






cc 


M(cc)[GeV] 


PS 


75 


'So 


0-+ 


Vc 


2.980(1) 


VC 


7m 


'S^ 


1"- 


J/^ 


3.097(1) 


SC 


1 


'Po 


0++ 


XcO 


3.415(1) 


AV 


757m 


'Pi 


1++ 


Xcl 


3.510(1) 



TABLE I. Charmonium states in different quantum number channels taken from the particle data book [30l] . 



The Euclidean temporal correlation function G(t, p) can then be defined as 

Gff(T,p)= J d^xe-'P-" (Jh(t,x)Jh(0,O)), 

where Gh(t, p) is the analytic continuation of D^(t,p) from real to imaginary time 
Gh(t,p) = Z?+(-zt,p). 



(11) 



(12) 



By using the KMS relation and the above equation, one can easily relate the correlation function to the spectral 
function, 



poo 

Jo ^'^ 

where the integration kernel K{uj, t) is 
cosh(u;(T - 1/2T)) 



sinh(a;/2T) 



(13) 



(14) 



Note that the kernel K{uj, r) is symmetric around t = 1/2T. 

Because of asymptotic freedom the spectral functions at very high energy are expected to be described well by the 
pro pag ation of a free quark antiquark pair. In this noninteracting limit the spectral function is analytically given 

by 



N, 



PH(^) = ^e(^^-4m^)a.^tanh — Wl- — 



AT J 



2m 



(1) (2) 



2m 

LU 



Nr. 



LU (S(w), 



with 



h=- 



d^k dnpiioi^) 



d'^k 97ii?(wk) 
2tt'^ duj\r 



27r2 ujI 



(15) 



(16) 
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FIG. 1. Left: free spectral functions on the isotropic lattice = 1) versus free spectral functions in the continuum limit. 
The vector spectral functions pvc{oj)/uj^ are plotted as function of ui/T for A''r = 24, 32, 48 at fixed T. Right: the free lattice 
spectral function on the isotropic lattice = 1) versus the free spectral function on the anisotropic lattice = 4). The spectral 
functions pvc{^)/ijj'^ are plotted as function of u/T for A'^t = 48 at fixed T. In both plots the value of quark mass m by T is 
fixed to be 3. 



The coefficients a^'^''^-' can be read off from Ref. ^ij, o;^ = + and np{iij]^) is the Fermi distribution function. 
Note in the above expression that there is a term proportional to uiS{uj), implying a r independent contribution to the 
correlation function. This contribution is also known as a zero mode contribution [s^ . For correlators with massive 
quarks, the zero mode contribution vanishes only in the PS channel. 

On lattices with finite temporal extent Nt the spectral functions suffer from lattice cutoff effects. As shown 
in the left plot of Fig. [TJ the free lattice spectral functions for the Wilson fermion discretization, which is used 
in this work on isotropic lattices of temporal extent Nr = 24, 32, 48, strongly deviate from the free continuum 
spectral function in the large oj region. In contrast to the continuum case the lattice spectral function starts from 
^min/T = 27VT-log(l + ma^/^). It has two cusps at 2Ar^log(l + (2 + ma„)/C) and 27V^log(l + (4 + ma„)/C) and 
vanishes at ujmax/T — 2Nr log(l + (6 + maa)/£,) [l7| . Here maa is the value of the quark mass m in units of spatial 
lattice spacing ao-, ^ is the anisotropic factor, i.e. the ratio of lattice spacing in the spatial direction over that in the 
temporal direction, ^ = aa-fur- In Fig. [T] the value of quark mass m/T = mNracr /£,=3. These lattice cutoff effects 
can be well separated from the region of physics interests when the number of points in the temporal direction Nr is 
large, i.e. the lattice spacing at fixed temperature T ~ l/Nrar becomes small. In this work we use isotropic lattices, 
i.e. ^ = 1. In order to increase Nr, which is very crucial in the spectral function analysis, an economic way is to 
perform simulations on anisotropic lattices. For instance, the anisotropic factor ratio ^ = 4 or larger has typically 
been used in previous calculations (TT| - [l3l [l6j . However, the lattice spectral functions on anisotropic lattices > 1) 
are much more distorted. As seen from the right plot of Fig. [U which shows the free lattice spectral function with 
Nr = 48, the lattice spectral function on the ^ = 4 anisotropic lattice vanishes at a much smaller energy, almost half 
of that on the isotropic lattice. As a consequence, the two cusps move closer to the region of physics interests. Thus 
lattice spectral functions obtained on anisotropic lattices are more contaminated by lattice cutoff effects, e.g. lattice 
simulations with Nr on the = 4 anisotropic lattices roughly correspond to those with Nr/2 on isotropic lattices. 
We therefore prefer to work on an isotropic lattice, although it is much more time consuming to generate gauge field 
configurations. 

At finite temperature, i.e. in the interacting case, due to the conservation of the vector current, the spectral function 
Pqq in the 7*^ channel will contribute a r independent constant to the correlator 

Poo = '^'^Xoo(^S{uj), (17) 
G^o = Txoo- (18) 

In the spectral function p]^ in the 7* channel on the other hand the ujd{uj) contribution present at infinite temperature 
changes into a smeared peak at finite temperature. From linear response theory the shape of this peak is expected to 
be a Breit-Wigner like distribution [sl] 



pJ;(w«T) = 2xoo|t^^, V = j^- (19) 
Ai Lo'^ + AID 
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Here M is the mass of the heavy quark, 77 is the drag coefficient and D is the heavy quark transport coefficient defined 
in Eq. ©. The contribution from Eq. ([T^ to the correlation function is generally called the smeared zero mode 
contribution. 



III. DETAILS OF LATTICE SIMULATIONS 



In this work we present results based on quenched lattice QCD simulations performed on isotropic lattices using 
C'(a)-improved Wilson (clover) fermions. The simulation parameters are shown in Table HH As we are interested 
in temporal correlation functions there is a need for very fine lattices in order to have enough data points in the 
temporal direction and reduce lattice cutoff effects. We thus performed simulations on lattices with lattice spacing 
ranging from O.OIfm to O.OSfm corresponding to the bare gauge couplings P — Q/ = 7.793, 7.457 and 6.872. At these 
(3 values the lattice spacing has been determined from the string tension parameterization with T^^j ^fa = 0.630(5) 
and y/a = 428 MeV [13, |3l|. The simulated temperatures range from about 0.75 Tc to 3 Tc- Simulations have 
been performed at T ss 0.75 Tc and T w 1.5 Tc with three different lattice spacings. This allows an estimate of the 
magnitude of lattice cutoff effects. For the higher temperatures, T k, 2.2 Tc and T w 2.9 Tc, simulations arc done only 
on the finest lattice. The number of correlator data points in the temporal direction is more than doubled compared 
to the finest lattice used in our previous study [l^. To reduce the volume dependence on such fine lattices, we use a 
large spatial lattice A^o-=128. The spatial extent thus ranges from 1.3 fm on the finest lattice to 3.9 fm on the coarsest 
lattice, which is in all cases significantly larger than the charmonium diameter. A subset of these lattices was used 
previously for the study of light meson spectral functions in Ref. (36j . 



P a[fm] 


a-i[GeV] L<,[fm] 


csw 






T/Tc ^ conf 


6.872 0.031 


6.43 


3.93 


1.412488 


0.13035 


128^ X 32 


0.74 


126 












128=* X 16 


1.49 


198 


7.457 0.015 


12.86 


1.96 


1.338927 


0.13179 


128^ X 64 


0.74 


179 












128^ X 32 


1.49 


250 


7.793 0.010 


18.97 


1.33 


1.310381 


0.13200 


128^ X 96 


0.73 


234 












128^ X 48 


1.46 


461 












128^ X 32 


2.20 


105 












128^ X 24 


2.93 


81 



TABLE II. Lattice parameters and number of configurations used in the analysis with a clover improved Wilson fermion action. 



P K He arrib T/Tc Nt aruAwi mRGi[GeV] m-^{m)[GeV] 



6.872 


0.13035 


0.13497 0.13130 


0.74 


32 


0.13305(2) 


1.592(4) 


1.255(2) 








1.49 


16 


0.13305(2) 


1.592(4) 


1.255(2) 


7.457 


0.13179 


0.13398 0.06201 


0.74 


64 


0.065430(6) 


1.4742(3) 


1.1739(2) 








1.49 


32 


0.065352(4) 


1.4734(8) 


1.1733(6) 


7.793 


0.13200 


0.13346 0.04143 


0.73 


96 


0.044245(7) 


1.358(3) 


1.093(2) 








1.46 


48 


0.044222(2) 


1.357(2) 


1.094(1) 








2.20 


32 


0.044280(6) 


1.359(3) 


1.096(2) 








2.93 


24 


0.04420(1) 


1.357(3) 


1.095(2) 



TABLE III. Quark masses on available lattices. Here rrif, stands for the bare quark mass, ttiawi is obtained from the Axial 
Ward Identity at the scale of /i = 1/a and rn^{m) denotes the renormalized quark mass in the MS scheme at the scale of 

All gauge field configurations were generated using a heat bath algorithm combined with 5 over-relaxation steps, 
whereby neighboring configurations are separated by 500 sweeps. For the fermion part the C(a) nonpcrturbatively 
improved Sheikholeslami-Wohlert action [37| has been implemented in our simulation with nonperturbativcly deter- 
mined clover coefficients, csw [ssj . listed in Table HIl The inversion of the Dirac matrix is carried out by using the 
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Conjugate Gradient (CG) algorithm. At the smallest lattice spacing, i.e. at ^ = 7.793, we measured two-point corre- 



128^ X 48, 128^ X 32 and 128^ x 24 corresponding to temperatures 0.73 T^, 
Since the temporal extent of the lattices is large and the exponentially 



lation functions on lattices of size 128'^ x 96, 
1.46 Tc, 2.20 Tc and 2.93 T^, respectively, 
decreasing correlation functions consequently become small at large distances, a rather stringent residue of 10^"^** in 
the CG algorithm has been implemented in our simulations. 

The nonperturbatively improved clover action used in our calculations removes 0{a) discretization errors. However, 
in calculations with heavy quarks, discretization errors of order am can also be large. We have estimated the quark 
mass values using the Axial Ward Identity (AWI) to compute the so called AWI current quark mass, mAwi, and the 
related Renormalization Group Invariant quark mass torgI; for the different lattice data sets [39l - l4l| . Here we used a 
nonperturbatively improved axial- vector current with coefficient ca taken from Ref. [38j . The commonly quoted quark 
mass for the heavy quark is the mass at its own scale. We therefore scaled ttirgi in the MS scheme to the scale fi = m, 
where the evolution of mj^{fi) to fi is done using perturbative renormalization group functions known with four-loop 
accuracy f42|-|44| . The resulting quark masses are listed in Table IIIII The mAwi is independent of temperature, which 
consequently makes also the RGI quark mass mRci and running quark mass ?Tij^(m) temperature independent. On 
the finest lattice am^j^^fi = m) is around 0.06 and thus the discretization errors proportional to the quark mass 
should be small. After adjusting the quark mass parameters for our calculations it a posterior turned out that the 
J/ip mass on the finest lattice is around 10% larger than the physical J/tJj mass. In the other two cases our choice of 
parameters reproduces the J/ip mass very well. 

As a local current, Eq. (|10p . is used in our calculations, it needs to be renormalized. 



j: 



H 



2nZHia,m,n = l/a)J'^'a-^ . (20) 
The renormalization factors Znia, m, /i = 1/a) are estimated using one- loop tadpole improved perturbation theory 
Zniamq^g^) = Zuiamq = Q,g^,a^ = 1) (l + bH{g^)amq) , (21) 



where Zniflmq = 0,f/^j^, a/.t = 1) are the renormalization constants for the massless quark case. Zniamq = 
O^gl—^afi = 1) has been determined perturbatively with two-loop accuracy for all the channels f45l - l47| and non- 
perturbatively for vector and axial- vector channels [i^. The coefficients ^Hif/jj-g) can be expanded in powers of the 
gauge coupling. 



bnig^) = 1 + Cphng^, 



MS' 



(22) 



These coefficients have been calculated at one-loop level [49|, |50| and in particular, bn for the vector channel has 
been determined nonperturbatively (48j . The resulting renormalization factors used in our calculations are given in 
Table HVl 



/3 


Zsc 


Zps 


Zvc 


Zav 


6.872 


0.92 


0.98 


0.97 


0.99 


7.457 


0.87 


0.93 


0.92 


0.93 


7.793 


0.87 


0.92 


0.91 


0.92 



TABLE IV. Renormalization constants of local operators for different channels. 



To check the magnitude of discretization errors, we analyzed the dispersion relation of the mesons. At nonzero 
"momentum" (pj^ ^ or a;„ ^ 0), the exponential drop of the spatial correlator may be described by an energy Esc 

^2 

G{z, p±,uJn) - exp{~Escz), El^ = + _| + m^^, (23) 

where w„ ~ 2'nnT arc the Matsubara frequencies, pj^ is the transverse momentum, and mgc is the screening mass 
which can differ from the pole mass if A{T) ^ 1. It is worth noting that the above ansatz is based on the dispersion 
relation in the continuum limit. 

We show the dispersion relation of the screening mass in the PS channel in Fig. [51 The results are obtained from 
calculations performed on 128^ x iV^ lattices at 0.74 {N^ = 64) and at 1.49 {Nr = 32) with a"^ = 12.86 GeV. 
Labels in the figure indicate whether spatial ("x") or temporal components ("r") of (pxjPyjPr) were chosen to be 
nonzero. The lines denote the dispersion relation obtained by fitting with an Ansatz of E^^ijj) = ap'^ + b. At 0.74 
Tc, for the results from the spatial directions, we have a good fit with parameters a=1.02±0.01 and &=9.530±0.013. 
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FIG. 2. The dispersion relation of the screening mass in the PS channel obtained hom lattices with /3 — 7.457. Labels in the 
figure indicate whether spatial ("x") or temporal components ("r") of (pxyPyyPr) were chosen to be nonzero. The lines denote 
the dispersion relation obtained by fitting with the form of Eq. (|23[) . 



Mass in GeV 

P J/ii Vc Xci xco 

6.872 3.1127(6) 3.048(2) 3.624(36) 3.540(25) 
7.457 3.147(1)(25) 3.082(2)(2f) 3.574(8) 3.486(4) 
7.793 3.472(2)(ff4) 3.341(2)(104) 4.02(2)(23) 4.52(2)(37) 



TABLE V. Meson masses (in GeV) for different charmonium states. The errors in the first bracket are statistical errors and 
the errors in the second bracket are systematic errors from effects of the physical distance. 



The applicability of the Ansatz ap^ + b indicates that our lattice is very close to the continuum limit; for the results 
from the temporal direction, even though we only have 3 data points, at this temperature, the data points have 
the same behavior as those from the spatial direction. We also performed a fit and obtained a=1.01±0.03 and 
6=9.539±0.033. The slope parameter a here equals A^^ in Eq. (|23p . The proximity of a to 1 confirms that at 0.74 Tc 
the screening mass is a good approximation for the pole mass. The meson masses obtained from the spatial correlation 
functions at T < Tc are shown in Table IVl When going to the higher temperature of 1.49 Tc, the data points from 
the temporal direction differ strongly from the fitting line for the results from spatial directions. Thus, the temporal 
direction is distinguished from the spatial direction and the breaking of Lorentz symmetry is clearly observed at this 
temperature. We also note that the screening mass at 1.49 Tc is about 10% larger than the mass determined in the 
confined phase. 

In the following sections, we restrict ourselves to the case of vanishing momentum and suppress the p indices. 



IV. EUCLIDEAN CORRELATORS ABOVE n 



Following Ref. [l3| we introduce the "reconstructed" correlator at temperature T from a spectral function determined 
at temperature T' 



POO 

Grec(T, T; T') - / duj K{t, T, uj) p{uj, T'). 
Jo 



(24) 



9 



Giec{T, T; T') at the temperature T is computed from a spectral function at the temperature T' and an integral kernel 
K at the temperature T. In the following subsections, we will study the Euclidean correlation functions at T > T^, 
with respect to a reconstructed correlation function that uses a spectral function obtained at T' < T^. 



A. Remarks on the reconstructed Euclidean correlation function 



We want to compare correlators calculated a± T > Tc to those at T < T^, i.e. we consider the ratio of the measured 
correlator to the reconstructed correlator [l^ , 

G,,,iT,T;T') du K{t,T,u) p{u:,T'y ^ ' 

This reduces the influence of the trivial temperature dependence of the kernel K{t^ T, uj) in the correlation function. If 
the ratio is equal to unity at all distances it would suggest that the spectral function does not vary with temperature. 
In fact, in order to obtain the reconstructed correlator at temperature T from a spectral function at T' one does not 
require any knowledge of the spectral function at that temperature. It suffices to know the correlator at T'. 

To arrive at the desired correlation function at temperature T, we first exploit the following relation [20| which is 
a generalization of the relation derived in Ref . [5l| , 



cosh[ujif-Nr/2)] ^ cosh[c.(f- - iV;/2)] 

sinh(a;iV./2) ,,^,.Xf'=N^ sinh(c^iV;/2) ' ^ > 

where T' = (a7V;)-\ T = (aTV^)"!, f = (r'/a) e [0, TV^ - 1], f = (r/a) € [0, iV^ - 1], iV;=m7V^, m e Z+ . Nr 
and N!^ are the number of time slices in the temporal direction at temperature T and T', respectively; f denotes the 
time slice of the correlation function at temperature T while f ' denotes the time slice of the correlation function at 
temperature T'. The sum over f ' on the right hand side of Eq. (|26p starts from f' = f with a step length of Af' = Nt 
and ends at the upper limit N'^ — Nr + f. After multiplying both sides of Eq. p6| with p{uj, T') and performing the 
integration over a;, one immediately arrives at 

G,,e(T,T;T')= J2 G{f\T'), (27) 

f'=f; Af'=Nr 

which shows that Grcc{T,T;T') is obtained directly by using the correlator G{t',T') at T' . Using relation ((27)) 
we can calculate G^cc{ttT]T') directly from the correlator data at temperature T'. An immediate consequence 
clearly is that one has a better control over systematic errors in the calculation of ratios used in Eq. In the 

following subsections, we will implement Eq. (j27p to calculate the reconstructed correlators and compare them with 
the measured correlation functions. We will discuss what can be learned about the modification of spectral functions 
from the analysis of correlation functions. In the following sections we will suppress the index T' in the Groc- 



B. Ratios of G{t,T) to Grec(r,r) 



We first investigate the temperature dependence of the pseudoscalar correlators. We show the numerical results 
for G{T,T)/Grec{T,T) at 1.46 T^, 2.20 and 2.93 on our finest lattice in the left plot of Fig. |3l Grec{T,T) are 
evaluated from the correlator data at T' = 0.73 Tc using Eq. (P7|) . Note that the error bars shown in the plots are 
statistical errors obtained from a Jackknife analysis. As seen from the left plot of Fig. [3l the ratio G(r, T)/Grcc(T, T) 
approaches unity at small distances and starts to deviate from unity at larger distances. Temperature effects start 
to set in at about 0.06 fm at 1.46 Tc and make the ratio smaller than unity. Deviations are about 5% at the largest 
distance. The small temperature dependence of the pseudoscalar correlator might indicate that the corresponding 
spectral function is subject to only small thermal modifications. When going to the higher temperature, 2.20 Tc, the 
temperature effects set in at a smaller distance (« 0.03 fm). The ratio rapidly drops and the deviation from unity 
(~ 8%) becomes larger at the largest distance r = 1/2T. Turning to the highest available temperature, i.e. 2.93 Tc, 
the temperature effects also set in at around 0.03 fm and the deviation of the ratio from unity increases to about 12% 
at the largest distance. This may suggest considerable modifications of the lowest state in the PS channel at this 
temperature. 
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FIG. 3. The ratio G{t,T) /Gicc{j,T) for PS (left) and Vu (right) channels as a function of the Euclidean distance r on our 
finest lattice with /? = 7.793 (a = 0.01 fm) at T = 1.46, 2.20 and 2.93 Tc. The reconstructed correlator dec is obtained directly 
from correlator data at 0.73 Tc- 



The ratio G(t, T) / Grcc{T, T) for the vector correlator Vu (summing over spatial components only) on our finest lat- 
tice is shown in the right plot of Fig. |3l Clearly this plot shows that the temperature dependence of G'(r, T) /Groc(T, T) 
is quite different from that in the PS channel. At all temperatures the ratios are larger than unity. This is already an 
indication that different temperature dependent contributions arise in the vector channel at large distances, related 
to the low frequency region in the spectral function. The temperature effects set in at larger distances compared to 
pseudoscalar correlators: around 0.1, 0.08 and 0.06 fm at 1.46, 2.20 and 2.93 respectively. A unique feature seen 
in the Vu channel is that the magnitude of G{T,T)/Grcc{T,T) at the largest distance r ~ 1/2T does not vary with 
temperature. All ratios deviate from unity by about 16%. In fact, the ratios seem to be to a good approximation 
a function of tT only. However, one has to be careful with the interpretation of this result in terms of bound state 
modifications as their effect may be compensated by possible positive diffusion contributions in the Vu channel at 
temperatures above Tc. We will examine this in more detail in the next two subsections. 
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FIG. 4. Same as Fig. [3] but for P wave states. The left plot is for the SC channel and the right one is for the An channel. 



The numerical results for the ratio G/Grcc for P wave states obtained on our finest lattice are shown in Fig. S) 
The left plot of Fig. |4] is for the scalar channel while the right plot is for the axial- vector channel. They show similar 
features as we have seen in S wave correlators: at short distance the ratio is close to unity while at the large distances 
the ratio deviates from unity and the deviations start at shorter distance at higher temperatures. We find a significant 
deviation of G/Groc from unity in both channels already at 1.46 T^: at the largest distance r = 1/2T, G/Grcc reaches 
about 1.9 in the SC channel and about 2.5 in the An channel. This deviation at the largest distance is much larger 
compared to the case of S wave correlators and the magnitude of this deviation at the largest distance decreases with 
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increasing temperatures. In order to connect these features with the thermal modification of bound states, one needs 
to separate the contribution of the smeared zero mode at low frequency, which is present in SC and An channels. 

We also compared the ratios G{t,T) /Grcc{T,T) obtained from the finest lattice with those from the two coarser 
lattices, ^ = 6.872 (a = 0.031 fm) and (3 = 7.457 (a = 0.015 fm). We found that the lattice cutoff effects are 
small in these ratios. The ratios at the largest distance obtained on the finest lattices are about 7% larger than 
those on the coarser lattices. This is mainly due to our choice of quark mass parameters, which lead to somewhat 
larger charmonia masses on our finest lattice. In fact, the exponential decrease of the ratio of correlation functions, 
G/Grcc ^ exp(— AMt), is controlled by the difference of effective meson masses M^s below and above Tc, i.e. 
AM ~ Meff(r > Tc) - Mcff{T < Tc). As the mass Meff(T ~ 0.75 Tc) is larger on our finest lattice, the high 
temperature mass M^tfiT ~ 1.5 Tc) differs less from the low temperature value and thus leads to a smaller value of 
AM. This in turn leads to a smaller decrease in G/Gjcc relative to the result on the coarser lattice and explains the 
somewhat larger values for G/Groc on our finest lattice. 



C. Smeared zero mode contributions 



As discussed in Section|TTl there are zero mode contributions in An, SC and Vu channels in the infinite temperature 
limit. In the study of spectral functions at temperatures below Tc, which will be shown in Section IVBl we found 
no zero mode contributions in SC, PS and Vu channels at 0.73 Tc- Thus smeared zero mode contributions (or 
smeared 6 functions) are expected to arise in An, SC and Vn channels only at temperatures above Tc- In An, SC 
and Vu channels, the information on bound states and smeared zero modes are thus entangled in the low frequency 
region and may partly compensate each other. In order to retrieve reliable information on possible bound states 
it is therefore necessary to filter out or separate the smeared zero mode contribution. To emphasize this point 
we note that previous studies suggested that the temperature dependence of G/Grec is in fact mainly due to zero 
mode contributions [3^ [52l . [53| . On the other hand, the smeared zero mode contribution is interesting in its own. 
For example in the vector channel it is related to the diffusion (process) of a single quark in the medium. In this 
subsection we will discuss the evidence we have for smeared zero mode contributions and thermal modifications of 
bound states at the correlation function level in more detail. 

To get a better understanding of the r dependence of the ratios defined in Eq. (|25p , we perform a Taylor expansion 
of the correlators at the largest distance accessible at finite temperature 



CJO 



G{r,T) =/ ^^pj.) co^M^(;- V;^)) (28) 
^ ' ' 'o 2tt ' sinh(w/2r) ^ ' 



2tt sinh(a;/2T) 
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This allows us to explore the properties of the low frequency behavior of the spectral function. Here we define the 
Taylor expansion coefficients, i.e. the time derivatives of the Euclidean correlation functions. 



1 d"G(T,T) 



dirTy 



tT=1/2 



1 dio -"H (29) 



27r \TJ sinh(w/2r)' 



as thermal moments [36j . By going to higher order thermal moments, one probes higher frequency region in the 
spectral function. In particular the value of the zeroth order thermal moment G^^-* is the same as the value of the 
correlator at the symmetry point, G{tT = 1/2). We have extracted the zeroth and second order thermal moments 
(G^^-* and G^^-*) from correlation functions on the finest lattices. The results are shown in Table IVTl 
We rewrite the Taylor expansion of the Euclidean correlators G(r, T) (Eq. (|28p ) as 

G(r, T) - G^") 5] i?2„,o f rT - - j , i?„,„ ^ ^ . (30) 

n— ^ 

Thus the ratio of measured correlator to the reconstrueted correlator can be expanded as 



which shows that the sign of R^'^ — iJ^oc determines whether G{t,T)/G-ccc{t^T) is decreasing or increasing with tT 



at large distances. Take the S wave states for example. From Table rVTl we find that R ' ' ~ R^^^ is negative in the 
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Vii channel and positive in the PS channel at all three temperatures above Tg. It indicates that G{t,T) / Giec{T,T) 
increases with tT in the Vn channel and deceases with tT in the PS channel at large distances at T > T^, which is 
consistent with Fig. [3l 



channel T/T^ G<°Vr=' G["^/r^ G'^Vt^ G?JJT'^ AG'-'^^ /T^ AG(^Vr=^ AG^/T^ 

Vii 1.46 0.955(5) 0.829(6) 46.39(4) 46.43(7) 0.126(8) -0.04(8) 

2.20 1.81(2) 1.561(9) 57.3(2) 59.2(1) 0.25(2) -1.9(2) 

2.93 2.33(2) 1.99(1) 59.6(3) 62.6(3) 0.34(2) -3.0(4) 

PS 1.46 0.858(8) 0.91(1) 44.73(3) 45.74(7) -0.05(1) -1.01(7) 

2.20 1.44(2) 1.56(2) 52.6(1) 54.2(1) -0.12(3) -1.6(1) 

2.93 1.68(2) 1.90(2) 51.2(1) 54.1(2) -0.22(3) -2.9(2) 

Au 1.46 0.708(6) 0.280(3) 23.80(4) 22.71(2) 0.428(7) 1.13(2) 45(4) 

2.20 1.57(3) 0.761(6) 39.9(2) 40.0(1) 0.81(3) -0.1(2) 

2.93 2.18(3) 1.186(7) 47.9(2) 48.8(3) 0.99(3) -0.9(3) 

SC 1.46 0.493(5) 0.259(5) 21.32(2) 20.15(3) 0.234(7) 1.12(1) 33(2) 
2.20 0.99(2) 0.665(7) 32.3(1) 32.21(9) 0.33(2) -0.1(1) 
2.93 1.26(2) 0.980(9) 35.2(1) 36.2(2) 0.28(2) -1.0(2) 



TABLE VI. Thermal moments extracted from correlator data on the finest lattice (/3 = 7.793). AG*-"^ is the diflference between 
G'"' and dec as defined in Eq. I[40\. 

As mentioned before the contribution from (smeared) zero modes and bound states to the correlation function is 
difficult to disentangle by investigating the ratios G{t,T) /G-[cc{t,T). To investigate the modification of bound states, 
one has to separate the zero mode contribution. Since the ijj5{lo) term only contributes to G*-"-*, one can construct 
quantities that only include the higher order thermal moments. One possibility is to look at the ratio of differences of 
the correlators at neighboring Euclidean time slices to the difference of the corresponding reconstructed correlators, 

G''^(T,r) ^ G(f,T)-G(f + i,r) 

Gdiff (r, T) - G,.ec(f , T) - G,ec(f + 1, T) ■ > 

This approximates the ratio of the time derivative of the measured correlators to the time derivative of the recon- 
structed correlators at f -|- 1/2. Here f = r/a and the difference of correlators G(r, T) = G(r, T) — G(f -I- 1,T) can 
be expanded as 

G(r,T) ^ G(i) ( - 2 + ^ ) ' ^ (33) 

n— ^ ^ ^ 

where 



1 d"G(r,r) 



n\ d(TTY 

tT=1/2-1/2JVt 

Thus Eq. (|32p can be rewritten as 



2 Alo /a;\" , , sinh(a;/2A^^r) 

nil 2^ (t) sinh(:./2T) ' ^''^ 



G'^^^ir, T) _ G(i) ^ ^^3, _ ^3,,^ (,T-\ + — ]+...]. (35) 



Alternatively one can consider the ratio of midpoint subtracted correlators 

G-(., T) _ G(r, T) - G(iV./2, T) _ G^^ , _ (rT-'^\■^. (36) 



G-b(r,T) G,.ec(r,r)-Grec(iV,/2,r) G?l\ ^ \ 2 

As seen from Eq. and Eq. the zeroth order thermal moment G^^^ drops out in ratios G'^'*(t, T)/Gf^{T, T) 
and G^^^{T,T)/Gf^{T,T). Since a W(5(a;) term in the spectral function only contributes to the zeroth order thermal 
moment and, moreover, its contribution vanishes in the higher order moments G^"-^' and G'-"-^\ it is thus possible 
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to completely remove the zero mode contribution in the two ratios G'^'^{T,T)/Gf'^^{T,T) and G'''^{t,T)/G;'^'^{t,T). 
However, at finite temperature above T^,, the lliSIlu) term is likely to be smeared out as a Breit-Wigner like distribution 
(Eq. (fT9|)). This Breit-Wigner like distribution in the very low frequency region of the spectral function does not lead 
to a T independent constant, and it contributes to the thermal moments at all orders. Thus the smeared zero mode 
contributions cannot be completely removed from the above two ratios. However, these contributions, which are 
located only in the frequency region w < T in the spectral function, are suppressed at higher orders of the thermal 
moments due to the presence of a factor (^)" in Eq. and Eq. 



1.050 



1.025 



1.000 



0.975 



0.950 



0.925 



" diff sub 

1 .46 Tj, ^-m-^ ^^H- 

2.20 Tj, 

2.93 ^-.-^ 



x[fm: 



1.050 



1.025 



1.000 



0.975 



0.950 



it 



0.05 



0.1 



0.15 



0.2 



0.925 
0.25 



1 .46 
2.20 
2.93 



diff sub 



x[fm] 



0.05 



0.1 



0.15 



0.2 



0.25 



FIG. 5. The ratio C^'^ {t,T) /Gf^^ {t,T) (G^''^{t,T)/G^,^^{t,T)) of S wave sates as a function of the Euclidean distance r on 
our finest lattice with ^ = 7.793 at T = 1.46, 2.20 and 2.93 Tc. The superscript "x" denotes eitlier "difF' or "sub". The left 
plot is for the Vu channel and the right one is for the PS channel. 



In Fig. Owe show results for G'^'^/G™^' and C^'^/Gfl^ in the V^^ (left) and also in the PS channel (right). The 
open symbols denote the ratio G^'^/G™^^ while filled symbols label the ratio C^'^/Gfi^. The ratios G^'VG?^,'' and 
G'^'^/Gf^^ give similar results at all distances. In the Vu channel we observe that values of G'^'^/G^^^ and G'^'^/G^^f 
are much smaller than those of G/G^^c large distances. The values at the largest distance are reduced by almost 
15%. In the PS' channels deviations of G^^^ /G^^ and G'^^^ /G'^^ from unity arc also reduced compared to those 
of G/Gicc at large distances. However, the change is not as large as in the Vu channel. The values at the largest 
distance are increased only by about 3% at both 1.46 and 2.20 Tc and about 6% at 2.93 Tc- The larger changes that 
occur in the Vu channel correlators relative to those in the PS channel may be understood in terms of large smeared 
zero mode contributions that contribute in the vector channel and get almost completely eliminated in subtracted 
correlation functions. Such contributions do not seem to be present in the PS correlator and the resulting changes 
are thus smaller. We will confirm this interpretation through the explicit construction of the spectral functions in 
the next section. We also note that in the subtracted correlators there are significant differences in the pseudoscalar 
and vector channels. In fact, in the subtracted correlators the situation now seems to be reversed compared to the 
unsubtracted correlators. At T = 1.46Tc the subtracted vector correlator now stays close to unity at all distances rT, 
while we observe a clear drop in the pseudoscalar correlator. Of course, one has to keep in mind that the subtracted 
correlators also modify the spectral contributions in the bound state region and thus may also suppress contributions 
that result from modifications of the bound states. The weak temperature dependence seen in Fig. [5] thus does not 
necessarily mean that the bound states suffer negligible thermal modifications. We will look into this issue more 
closely in the next subsection. At the two higher temperatures, T ~ 2.20Tc and 2.93rc, the subtracted correlators in 
both channels show almost identical behavior. 

In addition to the S wave states we also examined G'^'^/G™^'' and G'^'^ /G'^^ for the P wave states. The corre- 
sponding results are shown in Fig. [SJ where the left plot is for the Au channel while the right plot is for the SC 
channel. The magnitudes of the ratios for both Au and SG channels are greatly reduced compared to the ratios 
shown in Fig. |4l This behavior is quite similar to the ratios in the Vu channel and it suggests that the strong rise 
seen in G/Grec m Fig. 2] could be partly due to smeared zero mode contributions. However, as mentioned in the case 
of the vector correlator, also in the Au and SG channels the bound state contributions get modified in subtracted 
correlators. We will examine this in more detail in the next subsection. 
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FIG. 6. Same as Fig. [5] but for P wave states. The left plot is for the An channel and the right one is for the SC channel. 

D. Difference between G(r,T) and Grcc{r,T) 



The behavior of the ratios G{t,T) / G,c^{t,T), G^'^/G™^' and C^'^/C^^^ provides some insight into the relative 
importance of different frequency regions for the structure of correlators. However, the contribution from smeared zero 
modes and bound states to the correlation function is still difficult to disentangle at this point. Further information 
is gained by looking into the differences between measured correlators and reconstructed correlators 



AG(T,T)/r= = (G(t,T) - G,«(T,T))/r= 



j ^ Ap{u,.T)K(u,,rT)/T' 



where 



and 



ag(") = g(")-g1"). 



(37) 
(38) 

(39) 
(40) 



The difference between the measured correlator and the reconstructed correlator provides information on the difference 
of spectral functions below and above Tc- As can be seen from Eq. dnj) and Eq. the intercept (AG^")) and 

curvature of G{t,T) — Groc(T, T) at large distances (AG^^^) are related to the values of G(t, r)/Gicc(T, T) and 
G'^^^{t,T)/G^^^{t,T) at the largest distance, respectively. The values of AG'°^ and AG*^^', obtained by performing 
a two-parameter quadratic fit to (G(r, T) — Grcc{T,T))/T^ , are hsted in Table IVll 

We first show the differences {G{t,T) — G^cc{t,T)) /T^ for S wave states in Fig. [71 The first thing to notice is 
the change in the dependence of (G(t, T) — Grcc{T,T))/T^ on Euchdean time tT as one raises the temperature. In 
Vii and PS channels we observe that {G{t,T) — Grcc{T,T))/T^ increases with tT at all temperatures. This is also 
reflected in the differences of second order thermal moments AG^^-*, which are clearly negative at the two highest 
temperatures in the Vu channel and at all temperatures in the PS channel as seen from Table IVll It is obvious that 
AG(t, T) would decrease with tT and AG*^^-' thus would be positive, if Ap(a;) > for all w. Increasing differences 
of correlators AG(r, T) and negative values for AG'-^^ thus indicate that Ap{uj) is negative in some energy range. 
On the other hand, one clearly cannot draw the reverse conclusion, i.e. we cannot rule out that A/9(aj) is negative 
in some energy range even if AG(r, T) decreases with tT and AG^^-* is positive. However, in that case regions with 
Ap(a;) < need to be compensated by regions with an enhancement in Ap{uj) > 0. To this extent it is worthwhile to 
note that the presence of smeared zero mode contributions above Tc, which do not have a counterpart below Tc, will 
give positive contributions to Ap{uj), while disappearing bound states will lead to negative contributions to Ap(w). 
We thus conclude that there must be some energy regions in which Ap is negative in the Va channel at the two 
highest temperatures and in the PS channel at all temperatures we examined. At 1.46 the second thermal moment 
AG^^-' is slightly smaller than zero in the Vu channel. It is manifested in the behavior that G{t,T) — Grec(T, T) is 
almost flat at large distance and slightly increasing with distance at 1.46 Tc- However, the increase with tT is not 
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statistically significant. If modifications of Ap(a;) would arise from the smeared zero mode only, i.e. in the w < T 
region, its contribution to correlation function would be either a constant at all distances or decreasing with distances. 
The former case corresponds to a LdS{uj) term in Ap(w) and the latter case corresponds to a smeared Breit-Wigner 
like distribution. Fig. [7] thus indicates that some modifications in the u > T region of the Vu channel are likely to 
happen already at 1.46 Tc- Combining the above discussion with the information on smeared zero mode contributions 
gained from the analysis of the subtracted correlators discussed in the previous subsection (and the following section 
on spectral functions) we conclude that Fig. [7] provides strong evidence for modifications of the spectral functions 
of S wave states that lead to Ap(a;) < in some energy range for all T > lAGTc- Also note that the differences 
{G{t,T) — Grcc{T,T))/T^ are negative at all temperatures in the PS channel. It is in contrast to the positive values 
of {G{t,T) — G,-cc{t,T))/T^ in the Va channel. This too suggests the existence of a significant smeared zero mode 
contribution in the Va channel and larger thermal modifications in the bound states in the PS channel at T > T^- 
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FIG. 7. (G(t, T) - Grec(T, T)) /T^ of S wave sates as a function of the Euclidean distance r on our finest lattice with /3 = 7.793 
at T = 1.46, 2.20 and 2.93 T^. The left plot is for the Vu channel and the right one is for the PS channel. 



The results for the differences {G{t,T) — G^cc{t,T)) /T^ in An and SC channels are given in Fig. [51 At the two 
highest temperatures (G'(t, T) — G^cc{t, T))/T^ increases with tT in both channels. This indicates that Ap < in some 
energy region in these two channels at the two highest temperatures. However, the interpretation of the decreasing 
differences of correlation functions in and SC channels, which we observe at 1.46Tc, is a bit more complex. These 
correlation functions receive positive contributions to Ap{u;) from smeared zero modes and positive values for Ap^u) 
thus will arise at small w/T, i.e. for lo ^T. However, if in addition Ap{uj) would not change or stay positive also in the 
region w > T, the higher order moments AG^") would still fulfiU an inequality, AG^"' > 2AG(2) > 24AG''^' > 0. This 
inequality, however, does not hold in An and SC channels at 1.46Tc as is evident from Table IVTl We thus conclude 
that also in the P wave spectral functions at 1.46 Tc thermal modifications in the ui > T region occur that lead to 
negative Ap{uj). It is thus plausible that the P wave states disappear already at T = lA6Tc (as expected) but smeared 
zero mode contributions in these channels are so large that they still give the dominant contribution to the shape of 
AG(t, T) at this temperature. It is also worthwhile to note that the magnitude of difference, G(t, T) — Grcc{T,T), is 
smaller at 2.93 Tc than at 2.20 Tc in the SC channel. This too may be due to a partial cancellation of effects arising 
from the smeared zero mode contribution and those originating from modifications of the bound states. 



E. Charm quark diffusion coefficient estimated from correlations functions 

The charm diffusion coefficient is related to the smeared zero mode contribution in the Vu channel. As the very 
low frequency structure of the spectral function should manifest itself most strongly at the largest distance of the 
correlation function, the symmetry point of the correlation function G{tT = 1/2) should be strongly infiuenced by 
the transport contributions. 

At 1.46 Tc, one may assume that the intermediate and high frequency region of the spectral function is similar to 
that at 0.73 Tc- Also based on the fact that there is no zero mode contribution at 0.73 Tc in the Vu channel, one could 
then estimate the charm diffusion coefficient by fitting the value of G(l/2) — Grec(l/2). Here we ignore the difference 
in the intermediate and high frequency region at 0.73 and 1.46 Tc and only use the ansatz given in Eq. (|19|) for the 
transport peak to fit the value of G(l/2) — Groc(l/2) at 1.46 Tc ■ There is only one parameter, i.e. the heavy quark 
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FIG. 8. Same as Fig. [7] but for P wave states. The left plot is for the An channel and the right one is for the SC channel. 



mass M, that needs to be fixed to obtain the charm diffusion coefficient D. We note that the correlation function 
calculated from Eq. (jl9p at tT ~ 1/2 decreases faster with increasing heavy quark mass than with decreasing D. 
Thus there exist a maximum value of quark mass beyond which no solution for D exists. The maximum value of 
quark mass here is around 1.8 GeV. As the quark mass extracted from correlation functions is around 1.0 GeV (see 
Table HlH) we vary the charm quark mass from 1.0 GeV to 1.8 GeV. The charm diffusion coefficient D multiplied by 
2ttT then ranges from 0.6 to 3.6, i.e. 

M= 1.0 GeV, 27rrL>w0.6, (41) 
M = 1.8 GeV, 2nTD « 3.6. (42) 

If there is no negative contribution from A/9(w) to G{l/2) — G'roc(l/2), 3.6 could be an upper bound on 2ttTD at 1.46 
Tc- We also performed a fit with a linear form of btu describing the very low frequency behavior of the vector spectral 
function. Fitting to the difference of correlators at the symmetry point G{l/2) — Grcc(l/2) at 1.46 Tc gives: 

ttT 

2ttTD = 6 « 2. (43) 

3xoo 

Clearly the estimate of the charm diffusion coefficient is sensitive to the ansatz used for the fits. However, the charm 
diffusion coefficient estimated from these two different Ansatze are compatible. 

When going to higher temperatures at 2.20 and 2.93 Tc, the interplay between the change of bound states and 
diffusion part in the spectral function becomes complicated, thus it is not convincing that one may get a reasonable 
estimate of the charm diffusion coefficient by using a simple ansatz consisting of only a transport peak. Nevertheless 
in the following we will use these current estimates at 1.46 Tc as input for the default models that have to be supplied 
to the MEM analysis. Further details on the choice of default models are given in Appendix [X\ 



V. SPECTRAL FUNCTIONS 



In the previous subsections we found that the flatness and the small deviation from unity of the ratios (G/Groc(''", T), 
G'^'^/Gf^^{T,T) and G""""^ /G'','^^{t,T)) do not necessarily mean that thermal modifications of the ground states are 
negligible. The more relevant quantities to look at are the sign of AG(2) and the relative strength of AG^''^ , AG'^) and 
AG'*^^. However, from all these quantities, only a qualitative understanding of the thermal modifications of spectral 
functions at temperatures from below to above Tc can be deduced. To really explore the properties of charmonium 
states at different temperatures, one has to advance to a direct analysis of spectral functions by using the Maximum 
Entropy Method. As in the vector channel, the contributions from the bound states and the diffusion part are 
entangled, it would also be helpful to provide information on the diffusion part of the spectral function estimated in 
Section llV El into the default model in MEM analyses. 

In this section we start with a brief introduction to the Maximum Entropy Method in subsection IV Al We then 
discuss charmonium spectral functions below and above Tc from the Maximum Entropy Method in subsection IV Bl 
The estimation on the value of the charm quark diffusion coefficient from vector spectral functions at T > Tc will be 
given in subsection IV CI 
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A. Maximum Entropy Method 



Inverting Eq. (|13p to extract the spectral function is a typical ill-posed problem. At finite temperature the inversion 
is more complicated than at T = 0, since the temporal extent is always restricted to the temperature interval, 
< r < 1/T. The spectral functions we want to obtain are continuous while the correlators are calculated at a finite 
set of Nr Euclidean time points which typically are C(10). An infinite number of solutions thus exists. The task then 
is to select the most likely solution which is consistent with additional constraints. Because of the positivity and the 
normalizability of the spectral function it can be interpreted as a probability function. The guiding principle for the 
selection thus can be the Bayesian statistical inference, which is the basis for the Maximum Entropy Method. 

The Maximum Entropy Method (MEM) is a widely used tool for extracting spectral functions from correlation 
functions. It was introduced to lattice QCD by Asakawa et al. [H^l and has been successfully applied to lattice QCD 
data at zero temperature to extract the parameters of the ground state and excited states of hadrons (55l - [60l | . The 
application to finite temperature lattice QCD has also been explored [lOl - fTsl . [TgI. [l8l - [2ll . [6ll l62j . Based on the Bayesian 
theorem, MEM provides a way to select a unique spectral function p(uj) and transfers the problem of specifying a 
parameterization of into the problem of specifying a likelihood function and a prior probability. The most 
probable spectral function p(ijj), given lattice data G and prior information H, can be obtained by maximizing the 
conditional probability 

P[p\GH] = exp(a5[p] - L[p]), (44) 
where L[p] is the standard likelihood function and the Shannon- Jaynes entropy S[p] is defined as 

p{uj) 



p{uj) - m{uj) - p{u})\og 



m{uj) 



(45) 



Here m{uj) is the default model which introduces the prior information on the spectral function p(a;) as the input, 
e.g. p{uj) is positive-definite; a is a real and positive parameter which controls the relative weight of the entropy S 
and the likelihood function L. The final spectral function is expressed as an integral over a: 

p{uj) ^ j dapaP[a\G] I JdaP[a\G], (46) 

where P[q!|G'] is the posterior probability of a given data G and pa is the most probable spectral function for a certain 
a. 

As pointed out in Ref. [11], the integral kernel K{t,uj) diverges at vanishing w, 
9T 1 

K{T,c.) = — + {—-r + Tr')co + 0[co]\ (47) 
cj oi 

In order to explore the low frequency behavior of spectral functions, it thus is of advantage to introduce a modified 
kernel that is free of this divergence and leads to a redefined spectral function. Since K{t, lu) has the following 
property, 

Nr-l 

K{uj,t)^ l/tanh(cj/2), (48) 

T = 

and lim,.,_^no tanh(a;/2) — 1, wc implemented in our analysis the following modified version of kernels and spectral 
functions [l9l. [gsI] 

K{t,uj) = t£inh{uj/2)K{T,uj), (49) 
p{uj) = coth{uj / 2) p{uj). (50) 

The modified kernel K{t, lj) cures the instability of MEM at w w and reproduces the behavior of the original kernel 
in the large uj region. 

Note that we have to specify the default model m{uj) to extract the spectral function from the correlator data. Thus 
choosing the default model (DM) is an essential part of the MEM analysis. Therefore all available prior information 
needs to be included in the default model, as it can strongly affect the output spectral function if the quality of 
the correlator data is not sufficient. It is natural to choose a default model which reproduces the behavior of the 
spectral function in the large uj region. Note that Eq. (|15p describes the propagation of a free quark antiquark pair 
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in the continuum limit. On the lattice, the high frequency part of the spectral function is strongly distorted due to 
lattice cutoff effects [l3, HH , as seen from Fig. [1] in Section |IT1 Rather than growing as the free lattice spectral 
function vanishes above a maximal frequency. Thus when extracting the spectral function from the correlation function 
calculated on the lattice, it is reasonable to use the free lattice spectral function as the prior information in the MEM 
analyses. The Breit-Wigner distribution, Eq. (fT9)) . replaces the ijj6{uj) term and is also added into the default model 
for the very low frequency region in our MEM analyses at finite temperatures. 

B. Spectral functions at finite temperature 

In this subsection we will discuss spectral functions obtained from MEM analyses at temperatures below and above 
Tc- We will mainly focus on the results from our finest lattice, i.e. a=0.01fm with /3 = 7.793. When we analyze 
correlation functions using MEM, we fix the number of points in the frequency space to N^j = 8000 and the step 
length aAa;=0.0005, i.e. we fix aujmax ~ 4 or LOmax ~ 76 GeV. If not mentioned otherwise we use the free lattice 
spectral functions with quark mass am « 0.06 as part of default models in our MEM analyses, whereby the value 
am 0.06 corresponds to the value of 7Ti^(to) listed in Table Hill 

The lattice spectral function is subject to lattice cutoff effects. These show up in the short distance behavior of 
correlation function and manifest themselves in the large energy behavior of spectral functions as shown in Fig. [1] 
To reduce lattice cutoff effects we omit some correlator data points at small distances, i.e. we use f = r/a = 
4, 5,6, • • • ,Nt/2 in the MEM analysis. In addition we need to take into account a default model modification of 
the large energy part of spectral functions that arises from perturbative corrections to the free field behavior. When 
using a free continuum spectral function as ansatz this is usually done by multiplying the large energy part, which is 
proportional to w^, with a suitably chosen constant. This cannot be done with our ansatz for the default models where 
we use at large energies the free lattice spectral functions that are cut off at some maximal energy uimax- Instead 
we rescale the free lattice spectral function in the default model DM such that the correlator Gdm{t,T), calculated 
from the default model, agrees with the lattice data at r/a = 4, i.e. we demand GdmIt/o = A,T) / G{t / a = 4,T)=1. 
To suppress the large w rise, in general we plot the spectral function p{io) divided by cj^ as a function of w. 
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FIG. 9. Default model dependences of the output spectral functions in Vu (left) and PS (right) channels at 0.73 Tc on the 
128'^ X 96 lattice. The plots in the right panels are blowups of the low frequency region of the left panels. "DM"s are the input 
default models while "spf's are the corresponding MEM outputs. 

We first look into the left plot in Fig. [SI i.e. the Vu channel below Tc- We test three different default models, 
"DM1" is a rcscalcd free lattice spectral function, "DM2" is a rcscalcd free lattice spectral function supplemented 
with a resonance peak located in the low frequency region and "DM3" is a rescaled free lattice spectral function with 
a transport peak described by e.g. Eq. ()19p in the very low frequency region. In the very high frequency region 
(a; > 35 GeV), as seen from the left panel of the left plot in Fig. [3J the MEM output just resembles the behavior of 
the input default models. In the low frequency region, as seen from the right panel of the left plot in Fig. |9l spectral 
functions obtained from MEM have a unique form that is different from that of the input default models, i.e. results 
are independent of the default model and thus might refiect stable features of the spectral function. Wc found that 
the default model dependence of the first peak is very weak. Though there are some variations in amplitudes of peaks. 
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the peak location of the first peak is always the same, at a; « 3.48 GeV, which in turn is close to the value of the 
screening mass obtained from the spatial correlator quoted in Table|Vl Thus this peak can be interpreted as the bound 
state peak of J/V'- It remains stable and robust in MEM analyses performed with quite different prior information. 
However, the width of this peak cannot be directly interpreted as the width of J /ip due to the limited statistics and 
small number of data points in the temporal direction. The second and third peak in Fig. |9] could be a mixture of 
higher excited states or MEM artifacts due to the finite lattice spacing and limited number of correlator point^. The 
output spectral function "spf3" , which is obtained from the default model "DM3" , has no transport peak although 
such peak is implemented in the default model. Wc thus conclude that at this temperature there is no (smeared) zero 
mode contribution in the vector channel, or in other words the charm diffusion coefhcient is compatible with zero at 
0.73 Tf.. We performed the same analysis in the PS channel (the right plot of Fig. and find that results in the 
PS channel are similar to those in the Va channel. We did not observe a (smeared) zero mode contribution in the 
PS channel as well. To analyze the modification of bound states in the spectral function with temperature the PS 
channel is therefore a good candidate. 
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FIG. 10. Same as Fig. [9] but for SC (left) and An (right) channels. 



Spectral functions of P wave states are shown in Fig. (TU] The default model dependence of the output spectral 
functions is small in both SC and An channels. We observed stable P wave ground states. As in the case of PS 
and Vii channels, no (smeared) zero mode contribution is observed in the SC channel at 0.73 T^. In the An channel 
we always found that some remnant of the input default model is present at w « for all three spectral functions, 
although the remnant is not obviously seen in the right plot of Fig. [101 Because of the noise level, it is however not 
clear whether this originates from the zero mode contribution or insufhcient quality of temporal correlator data. 

In the following we present the results for charmonium spectral functions above Tc- In Appendix |X] we give a 
detailed analysis of the default model dependence of our results and also quantify the influence of distance windows, 
the lattice cutoff etc. We found that the outputs from various default models are all compatible with those obtained 
by using the free lattice spectral functions with an added transport peak as default models. Consequently all MEM 
results shown in the following are obtained by using such default models. Since we found some excess of the spectral 
function in the low frequency region in the An channel, we use in that case the massless free lattice spectral function 
with a Breit-Wigner like peak at low frequencies as the default model. 

We now focus on the statistical error analysis of the spectral functions. This is done by using the Jackknife method. 
In the literature statistical errors are often given on the mean of p{u)) over a certain w region in the spectral function 
plot [l^ [l3, [l^ [l^ . Here we rather perform MEM analyses on Jackknife blocks and calculate the Jackknife error 
of the amplitude of each point in the spectral function. MEM cannot reproduce the correct width of the resonance, 
however, it gives a stable and reliable peak locations of the spectral functions. We thus also estimate the statistical 
errors of the peak location of the first peak of the spectral function at 0.73 Tc and 1.46 Tc (see Table IVlII . A signal 
for the dissociation of charmonium states then is a shift of the peak location and the relative broadening of the peak 
at different temperatures. 

We show the statistical significance of output spectral functions at > 2 GeV in PS (left) and Vn (right) channels 
in Fig. 111! The shaded areas are statistical uncertainties of amplitudes of output spectral functions from Jackknife 



^ We discuss this in more detail in connection with Fig. I16l in Appendix [XI 
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FIG. 11. Statistical uncertainties of output spectral functions in PS (left) and Vu (right) channels at all available temperatures. 
The shaded areas are statistical errors of amplitudes of output spectral functions from Jackknife analyses and the solid lines 
inside the shaded areas are mean values of spectral functions. The horizontal error bars at the first peaks of spectral functions 
at 0.73 Tc and 1.46 Tc stand for the statistical uncertainties of the peak location obtained from Jackknife analyses. 



analyses and the solid lines inside the shaded areas are mean values of spectral functions. From the left plot of Fig.fTT] 
it is apparent that at 0.73 T^. the spectral function in the PS channel has large uncertainties in the amplitude at the 
point which corresponds to the ground state peak location in the mean spectral function. However, even at the lower 
end of the error band, the amplitude is still larger than the peak ampHtudes at the higher temperatures within the 
errors. We also show the statistical uncertainties of the first peak location of the spectral function at 0.73 Tc and 
1.46 Tc as horizontal error bars in the left plot of Fig. 1111 Unlike the large uncertainties shown in the amplitude of 
the peak height, the peak location of the ground state peak at 0.73 Tc is well determined. A Jackknife analysis yields 
ruri^ = 3.31(4) GeV (see Table lVII[) . At 1.46 Tc this peak is shifted by about 0.8 GeV to around 4.1 GeV, as is seen from 
Table IVlIl At 2.23 Tc there is hardly a peak structure that can be identified within the statistical uncertainties. At 
2.93 Tc the spectral function flattens further. Thus this picture, together with the systematic uncertainties discussed 
in the Appendix, suggests that r]c is melted already at 1.46 Tc- 

In the right plot of Fig. [11] we focus on the resonance part of the spectral function in the Vu channel. One sees that 
the peak location of the spectral function at 0.73 Tc does not have an overlap with the peak location of the spectral 
function at 1.46 Tc and the amplitudes between these two differ a lot (see horizontal error bars and also values in 
Table rvTI|) . At both 2.20 Tc and 2.93 Tc there is hardly any peak structure. Together with the study of systematic 
uncertainties discussed in Appendix [XI this picture indicates that also J/ip is already dissociated at 1.46 Tc- 
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TABLE VII. The locations of the first peaks in different channels obtained from MEM. Errors are estimated from the Jackknife 
analyses. The numbers for peak locations are in units of GeV. 

The statistical errors on P wave spectral functions arc shown in Fig. [T^l Here the results for the SC channel are 
shown in the left plot. When going to temperatures above Tc, the structure of the ground state peak is basically gone 
and results in a rather flat spectral function. This signals the melting of Xco at T > lAGTc. 

The right plot of Fig. [12] shows the result for the An channel. As temperature increases from 0.73 Tc to 1.46 Tc, it 
becomes apparent that the location of the fist peak is shifted to the larger energy region. The bump seen at 1.46 Tc 
becomes much broader at 2.20 Tc and flattens at 2.93 Tc- An enhancement of the small energy part (2 GeV ^ w ^ 
4 GeV) in the spectral functions is also observed at the two highest temperatures. This originates from our choice of 
the default model in the channel where we did not introduce a quark mass threshold in the free lattice spectral 
function. The systematic uncertainties arising from the choice of the quark mass cutoff have been discussed in the 
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FIG. 12. Same as Fig. [TT] but for SC (left) and An (right) channels. 



Appendix in connection with Fig. 1191 The change of structures of spectral functions in the An channels suggests that 
Xci is dissociated already at lAQTc. 



C. Charm quark diffusion coefficient estimated from spectral functions 
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FIG. 13. Left: statistical uncertainties of transport peaks at T > Tc. Right: the resulting charm diffusion coefficients. The 
boxes stand for statistical error estimated from Jackknife method while the bars stand for systematic uncertainties from MEM 
analyses. The numbers for charm diffusion coefficients are listed in Table IVlIIl 



We now focus on the very low frequency part of the spectral function given in the vector channel, i.e. lo/T ^ 2.5 or 
w < 1 • • • 2 GeV at T/Tc ~ 1.5 • • • 3. The statistical uncertainties of the transport peaks observed in the vector channel 
are shown in the left plot of Fig.[T31 The statistical uncertainties on the amplitude of the peak are relatively small. The 
charm diffusion coefficient is related to the amplitude of the transport peak at vanishing frequency through the Kubo 
formula (Eq. ([9])). The current estimate for the charm diffusion coefficient D is summarized in the right plot of Fig.[T3l 
The boxes stand for the statistical uncertainties and the error bars reflect systematic uncertainties obtained from the 
analyses discussed in the Appendix. The bound for the systematic uncertainties for charm diffusion coefficients at 
all temperatures is obtained from the analysis of the default model dependence discussed in the Appendix and is 
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taken from the lowest and highest values in Fig. [Ml Fig. [151 Fig- [HI a-nd Fig. [191 The resulting numbers are listed in 
Table [Vini We find that the mean value of 2ttTD is around two at all three temperatures above T^. As the quark 
number susceptibility xoo increases faster with temperature than the amplitude of the transport peak, the mean 
values of 27tTD increases only slightly with temperature. The charm diffusion coefficient obtained at 1.46 Tc is the 
most reliable one among the three temperatures above Tc since more prior information is known at this temperature 
as discussed in Section llV El At 2.20 Tc and 2.93 Tc, due to the lack of precise prior information and a fewer number 
of data points that can be used in the MEM analyses, the uncertainties on the charm diffusion coefficient thus might 
be underestimated. 





27rTD 


xoo/r^ 


1.46 


1.8±0.7(stat.)tJ;^(sys.) 


0.20894(1) 


2.20 


2.0±0.4(stat.)t?-J(sys.) 


0.46900(2) 


2.93 


2.3±0.4(stat.)l?;?(sys.) 


0.66112(4) 



TABLE VIII. Charm diffusion coefficients and quark number susceptibilities xoo above Tc. The "stat." stands for statistical 
errors estimated from Jackknife method and the "sys." denotes systematic uncertainties obtained from MEM analyses in 
Appendix lAl 

To close this section, we have studied charmonium spectral functions in different channels at temperatures below 
and above Tc- The general properties of ground states at T < can be reproduced by the MEM analysis of temporal 
correlation functions. The extracted peak location of the ground state is close to the physical value and is very 
reliable. However, the width of the ground state peak cannot be reproduced in the MEM analysis with the current 
quality of the temporal correlator data. Thus the signature for the dissociation of charmonium states is the shift of 
the first peak location and relative broadening of the first peak. Comparing spectral functions below and above T^ 
our MEM analyses suggest that both the S wave states {rjc and J/ip) and P wave states (xco and xa) are melted 
already at 1.46 Tc- The charm diffusion coefficient is found to be compatible with zero at T < Tc and around I/ttT 
at our available temperatures above Tc. 

VI. CONCLUSIONS 

We have investigated the properties of charmonium states at finite temperature in quenched QCD on large isotropic 
lattices. The standard Wilson plaquette action for the gauge field and the nonperturbatively 0{a) improved clover 
fermion action for charm quarks were used in the simulation. In the current study lattices with three different lattice 
spacings were used to control cutoff effects in the charmonium correlators and spectral functions. Since the use of a 
temporal extent with a large number of Euclidean time slices is a very important ingredient in the current study, we 
calculated charmonium correlators on the finest lattices (a = O.Olfm) with relatively large lattice sizes of 128'^ x 96, 
128^ X 48, 128^ X 32 and 128^ x 24 at 0.73 Tc, 1.46 Tc, 2.20 Tc and 2.93 Tc, respectively. 

At T < 0.73 Tc we found stable and reliable ground state peaks of charmonium states from MEM analyses, in 
which the peak locations are almost the same as the corresponding hadron masses determined from the large distance 
behavior of spatial correlation functions at the same temperature. However, the width of the ground state peak still 
cannot be interpreted as the physical width of hadron states. Thus the dissociation of the ground states in the current 
study is signaled by the shift of the peak locations and the relative broadening of the width. At T > Tc we first 
calculated the reconstructed correlation function directly from the correlator data at T < Tc. The curvatures of the 
differences between the measured correlators and the reconstructed correlators indicate that there are obvious thermal 
modifications to spectral functions at T > 1.46 Tc in all channels. However, it is hardly possible to distinguish the 
zero mode contribution and the thermal modification of the bound states in the spectral function from the study at 
the correlator level alone. We then advanced to the analysis on the spectral functions using the MEM. We utilized 
an improved integral kernel to avoid the instability of MEM in the very low frequency region. We compared the 
output spectral functions from lattices with three different lattice spacings and concluded that cutoff effects are small. 
Results on our finest lattice, which are the most reliable ones, thus should not be affected by severe cutoff effects. 
Using the correlation functions on the finest lattices, we studied the variation of the output spectral functions using 
different default models both below and above Tc. We checked the systematic uncertainties arising from the number 
of data points used in the MEM analysis, the lattice cutoff effects present at short distances, and the dependence on 
the quark mass threshold of the free lattice spectral functions. Statistical errors of the spectral functions arc estimated 
using a Jackknife analysis. By comparing the spectral functions below and above Tc, our analyses suggest that both 
P wave states {xcO and Xci) and S wave states (J/ip and rjc) are dissociated already at 1.46 Tc. 
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The determination of dissociation of J/i/i and rjc already at 1.46 Tc is quite different from previous lattice QCD 
studies [Tol-[l6j. which predicted the IS charmonium states to be dissociated at T > 2 T^. Since most of the previous 
results are obtained on anisotropic lattices [TTI - [l6l | , lattice cutoff effects may strongly affect the physics deduced from 
the correlation functions. The variational method approach used in Ref. [3, [3 is a good way to enhance signals of 
hadron states which are well known to be there, e.g. the ground and excited hadron states at zero temperature, and 
will also contribute to a spectral decomposition of thermal correlation functions. However, the crucial question at 
finite temperature is not whether these states do contribute. Most important is the magnitude of such a contribution. 
This docs get modified in a variational approach and it is thus difficult to draw firm conclusions on charmonium 
melting from variational approaches. Several earlier studies based on the MEM analysis implemented an uniniproved 
integral kernel, which introduces an instability in the MEM algorithm in both low and high frequencies El, [3 ■ 
Most importantly, we have doubled the number of points in the temporal direction compared to our previous study 
on an isotropic lattice [loj . and have reduced the spatial lattice spacing by about a factor of 4 or more compared 
to studies performed on anisotropic lattices [ll| - [l6j . Comparing the number of points in the temporal direction in 
our study to the studies on anisotropic lattices, Nr in our study is about 1.5 times larger than that in Ref. [lll - [l5l | 
and compatible with that in Ref. [lal- Moreover, a very detailed MEM analysis on the default model dependences 
and the systematic/statistical uncertainties has been performed in the current study. All this supports a better 
control over systematic effects in our analysis that suggests that both S and P wave states disappear at T ^ 1.46 Tc- 
Obviously it does not necessarily mean that all the charmonium states are dissociated at the same temperature. The 
magnitude of thermal effects in the spectral functions is observed to vary in different channels. This may indicate 
that the charmonium states will dissociate at different temperatures. The fate of charmonium states at temperatures 
between 0.73 Tc and 1.46 Tc, however, remains unknown for us due to the lack of lattice data sets with appropriate 
temperature values in our study. A lattice QCD study of the screening masses extracted from spatial charmonium 
correlation functions suggests that rjc and J/ip may survive at T < 1.5 Tc in the hot mediu m |643. [65l . Furthermore, 
several lattice QCD studies of charmonia suggest that XcO and Xci melt just above Tc [I^, Il2l. Il3l|. The sequential 
suppression scenario [66l - l68j thus is not in contradiction to our results and, in fact, appears to be in accordance with 
the disappearance of bound states in the bottomonium "spectral function" from the latest experiment results [6^ . 
Thus lattice calculations of the temporal correlation function of charmonia at lower temperatures would be interesting 
and crucial to locate the dissociation temperatures. 

For the first time the charm diffusion coefficient has been estimated on the lattice directly from an analysis of 
spectral functions. The charm diffusion coefficient D is found to be compatible with zero at T < Tc and is about two 
times larger than the results from AdS/CFT calculations, l/27rT at T > Tc- However, more efforts are needed to 
reduce the current uncertainties on the charm diffusion coefficient obtained from lattice QCD calculations especially 
at higher temperatures. The pQCD results seem to approach to our findings when higher order corrections are 
included [25j . The heavy quark diffusion coefficients obtained from a T-Matrix approach with the internal energy, 
on the other hand, are close to our results [t^]- Recently another approach to calculate the heavy quark diffusion 
on the lattice [7l| has been suggested and an exploratory study has been carried out in Ref. [72[ . Further progress 
along this line has been made in Ref. (tsI W^ . in which the heavy quark diffusion constant multiplied by 2nT, i.e. 
2ttTD, is obtained to be in the range of (3. 5. ..5) at T « 1.5 Tc- It is quite impressive to see that this estimate 
for the heavy quark diffusion coefficient is close to our estimate of the charm diffusion coefficient, although totally 
different approaches have been used. Besides the Maximum Entropy Method, a Fourier method especially designed 
for addressing the low fre quen cy behavior of spectral functions has been introduced and has been used to estimate 
the electrical conductivity [75|, |76[. It would be interesting to implement this method also for the heavy quark sector. 

In this study the effects of dynamical quarks are not included. The general picture concerning the properties of 
charmonium states might not change significantly as concluded from the study of charmonium states in two flavor 
QCD [l^ . However, in the medium with sea quarks there exists a DD threshold and one might expect charmonium 
states to dissociate at lower temperatures. Moreover, it is difficult to predict how the dissociation temperatures are 
influenced by the change of the pseudo critical temperature T^ as Tc becomes smaller when dynamical quarks are 
included in the system [77|. For the charm diffusion coefficient, one knows that at sufficiently high temperature where 



pcrturbative QCD calculations are applicable 2ttTD becomes smaller when sea quarks are included [25[ . So far there 



are no studies that explore the charm diffusion coefficient in lattice QCD with dynamical quarks included. 
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Appendix A: Uncertainties of the spectral functions 

In the appendix we show the defauh model dependences and various systematic uncertainties of the charmonium 
spectral functions. For illustration we will only focus on the Va channel as the uncertainties in the other channels 
are very similar. The basic settings of MEM used here are the same as what we mentioned in Section |V] if without 
additional description. As mentioned in Section |Vl MEM cannot reproduce the correct width of the resonance but 
gives a stable and reliable ground state peak location of the spectral functions, so the signal for the dissociation of 
charmonium states is the shift of the peak location and relative broadening of the peak at different temperatures. 

The appendix is organized as follows: we will first show the default model dependences of output spectral functions 
in Appendix lA 11 The main point to check is that how the first peak location of the output spectral function changes 
when resonance peaks with different peak locations are provided in the default model. Then we will study the 
systematic uncertainties of the output spectral function from MEM in Appendix IX2l 

1. Default model dependences 

Because of the limited number of correlator data points in the temporal direction, the MEM analysis becomes more 
difficult at temperatures above Tc- As has been done at temperatures below Tc, to check the reliability of the output 
spectral functions from the MEM analysis, the default model dependence test is always the first thing one needs to 
do. In principle one should put as much physical information into default models as possible. This rule leads to a 
very straightforward default model dependence test for the spectral functions above Tc- That is to fully benefit from 
the two limits which we already know quite well: the free lattice spectral function at very high temperature and the 
spectral function obtained from MEM at a temperature below Tc- To put these pieces of information into the default 
model, one might be able to check to which limit, free or confinement limit, the output spectral function is closer. 
However, due to the fact that the spectral function at T < has a sharp ground state peak and the quality of 
correlator data at temperature above Tc is not sufficient, the MEM output basically reproduces the spectral function 
below Tc with negligible changes at all three available temperatures above Tc- Thus in the following default model 
dependence test we will not use the full information of the spectral function at T < Tc into the default model but 
rather the information of the peak location of the ground state ("DM2"s in the following analyses). Besides this, we 
will also use the free lattice spectral functions with some additional resonance peaks and/or transport peaks. The 
signature for the dissociation of resonances is then characterized by the shift of resonance peak location and relative 
broadening of the peak at different temperatures. The default modes used in the current study are summarized in 
Table im 



— — — channel 
default model — — — 




DM1 


flspf + BW 


DM2 


flspf + BW + resl 


DM3 


flspf + BW + res2 


DM4 


flspf + BW + res3 



TABLE IX. The default models investigated in the current section, "flspf stands for the free lattice spectral function, "res" 
stands for the resonance peak structure according to the relativistic Breit-Wigner distribution, "resl", "res2" and "resS" are 
resonance peaks with the peak locations equal, smaller and larger than the corresponding resonance peak locations in each 
channel at T < Tc. "BW" is a Breit-Wigner like distribution according to Eq. (|19|) and may vary in the width and amplitude 
in different default models. 

Because of the existence of a diffusion contribution to the Vu channel one has to check the dependences both on 
variations of the resonance part and diffusion part in the default models. Then for the default model dependence 
test we first fix the large u behavior of the default model by using the rescaled free lattice spectral function and 
vary the information on the very small uj part. i.e. the transport peak described in Eq. (|19p . The default models 
used here correspond to "DM1" ,"DM2" and "DM3" in Table Hxl without resonance parts. We show the result in 
Fig- [HI The upper panel of Fig. [T4l shows p{u>,T)/uj'^ as a function of uj in the large lj region while the lower panel 
of Fig. [HI shows p{uj,T)/{u!T) as a function of uj/T in the very low frequency region. "DM"s are the input default 
models while "spf's are the corresponding MEM outputs. In particular, the transport parts of "DM1" in the MEM 
analysis at 1.46 Tc are parameterized in the Breit-Wigner form with 2'kTD — 3.6 and M = 1.8 GeV obtained from 
the fit to the difference between the measured correlator and the reconstructed correlator by using a Breit-Wigner 
ansatz discussed in Section IIV El Because of the interplay between the contributions from the diffusion and resonance 
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FIG. 14. Default model dependences (varying the transport peak) of spectral functions in the Vu channel at temperatures above 
Tc. At each temperature the very large uj part of the default model is fixed to the behavior of the free lattice spectral function. 
Upper panel: p{uj,T) /uj^ as a function of lj, Lower Panel: a blowup of plots in the upper panel in the very low frequency 
region and plotted as p{u} ,T) / {loT) versus ui/T. "DM"s are the input default models while "spf s are the corresponding MEM 
outputs. 



parts it is difficult to make an estimate of the charm diffusion constant directly on the correlator level at the two 
highest temperatures. Here we simply apply the same value of charm diffusion D estimated at 1.46 to the default 
models at 2.20 and 2.93 T^- Looking at the output spectral functions at each temperature, we find that the variation 
of the very small uj part of the default model gives negligible effects to the intermediate uj part (resonance part) of the 
output spectral functions. Concerning the temperature dependence of the resonance peak, the upper panel of Fig. 1141 
shows that, already at 1.46 Tc, the ground state peak becomes much broader and its peak location is shifted to larger 
energies compared to that at 0.73 Tc (see Fig. [9]). When going to the higher temperature of 2.20 Tc one can hardly see 
a bump in the interesting oj region. At our highest temperature available, 2.93 Tc, we find that the large w part more 
or less resembles the shape of free lattice spectral functions and no peak structure is observed. For the transport peak 
shown in the lower panel of Fig. [T31 the prior information of charm diffusion D estimated from G(l/2) — (7^0(1/2) 
is put into very low frequency part of the default models. MEM shows the sensitivity to the very low frequency part 
and the output spectral functions differ from the default model in this very small energy region. We observe that the 
output transport peak has a weak dependence on the input default models at all three temperatures above Tc- And 
the amplitude of the transport peak at vanishing energy increases with temperature. 

After studying default model dependences by varying the transport peak in the default model on the output spectral 
function in the intermediate lo (resonance peak) region in Fig. 1141 we now fix the very low frequency (transport peak) 
part of the default model and vary the intermediate w (resonance part) behavior of the default models. The default 
models in the very low frequency part are fixed to have the same behavior as "DMl" in Fig.[T3]at each temperature. 
Again note that the transport part of "DMl" is parametrized as 2'kTD = 3.6 and M = 1.8 GeV as discussed in 
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FIG. 15. Default model dependences (varying the resonance part) of spectral functions in the Vu channel at temperatures 
above Tc- The transport part of the default model is fixed in each temperature same as that of "DMl"s in Fig. 1141 Upper 
panel: pluj ,T) / uj^ as a function of oj, Lower Panel: a blowup of plots in the upper panel in the very low frequency region but 
plotted as p{u ,T) / {uT) versus u/T. "DM"s are the input default models while "spf's are the corresponding MEM outputs. 



Section flV El We test four different default models as listed in Table HXl "DM1" is a rescaled free spectral function 
with a transport peak, "DM2" is a rcscalcd free spectral function with a transport peak supplemented with a resonance 
peak whose peak location is the same as that of the spectral function at 0.73 "DM3" and "DM4" are basically 
the same as "DM2" but with a resonance peak whose peak location is smaller and larger than that of the spectral 
function at T < Tc, respectively. We show the default models and their corresponding output spectral functions 
("spf's) divided by as functions of w in the upper panel of Fig. [151 At 1.46 Tc there is a minor default model 
dependence of the output spectral functions, but the trend is similar: the peak location is shifted to a location larger 
than the peak location of the spectral function at 0.73 Tc (peak location shown in "DM2") and the width becomes 
larger. At 2.20 Tc the default model dependence is a little stronger. This may be due to the smaller number of data 
points in the temporal direction and lower statistics. However, outputs from MEM still have unique differences from 
input default models and they all have a trend to resemble the shape of the free spectral function. At 2.93 Tc we 
have only 9 points in the analysis and together with the issue of the transport peak, the default model dependence is 
considerably stronger than that in the analysis at the other temperatures. Based on the results from 1.46 and 2.20 
Tc we do not expect the peak location of the resonance peak at 2.93 Tc shifts to smaller energies compared to the 
case at 0.73 Tc and would rather expect that the spectral function at this temperature is much closer to the spectral 
function in the noninteracting case. In the lower panel of Fig. [15] we enlarge the very low frequency part of the upper 
panel in Fig. [15] and show p(a;)/(wT) as a function of w/T. Unlike the case in the lower panel of Fig. [141 the change 
of the default model in the intermediate w part (resonance part) has a relatively large effect on the output in the very 
low frequency region. It could be mainly due to the compensation of the very low frequency part to the changes of 
corresponding resonance parts. Without the quantitative description of the transport peak wc can observe a trend 
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that the amplitude of the transport peak becomes larger with increasing temperature. 

In a short summary, in this subsection we have checked the reliability of the output spectral function in the Vu 
channel from the MEM analysis by varying the resonance part and the transport part of input default models. At 
1.46 Tc, the default model dependence of the resonance part is relatively weak and the resonance peak observed at 
T < Tc generally shifts to high frequency region and becomes much broader. At higher temperatures, the default model 
dependences of the resonance part becomes stronger due to the insufficient qualify of the temporal data, however, the 
general trend is that no clear peak structures are found and the spectral function gets closer to the noninteracting case. 
In the very low frequency region, we supplied the transport peak parameterized by the estimation in Section IIVEI 
MEM showed certain sensitivity to the transport peak. The transport peak in the Vu channel has weak dependence 
in the default models when only the transport peak in the default model is changed while it has relatively large 
default model dependence when only the resonance part in the default model is changed. The general trend is that 
the amplitude of the transport peak is increasing with the increasing temperature. 

2. Systematic uncertainties 

In this subsection we explore the systematic uncertainties of the spectral function from the MEM analyses. This 
study includes lattice spacing dependencies, a comparison of the spectral function below and above Tc for a same 
number of data points used to extract the spectral function, lattice cutoff effects at small distances and the dependence 
on the threshold of the free spectral function used in the default model. 
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FIG. 16. Output spectral functions from MEM in the Vu channel at temperatures below Tc from three different lattices: 
128^ X 32 with /J = 6.872 (a^^ = 6.43 GeV) at 0.74 Tc, 128=* x 64 with /3 = 7.457 (a^^ = 12.86 GeV) at 0.74 Tc and 128=* x 96 
with /3 — 7.793 (a~^ — 18.97 GeV) at 0.73 Tc. The small plot inside is the output spectral function in the whole energy region. 



First we look into the lattice spacing dependence of the output spectral function on our available lattices. We 
show spectral functions from the Vu channel at temperatures below Tc in Fig. 1161 The results are obtained from the 
lattices with a'^ = 18.97 GeV (/3 = 7.793, 128^ x 96), a"! = 12.86 GeV (/3 = 7.457, 128^ x 64) and a'^ = 6.43 GeV 
(/3 = 6.872, 128'^ x 32). The plot shows the behavior of spectral functions in the low frequency region (2 < w < 10 GeV) 
while the small plot inside shows the behavior in the whole frequency region. One can observe that with smaller lattice 
spacing the lattice cutoff effects (the cusps in the high frequency region) can be well separated from the physically 
interesting frequency region. As seen from the low frequency region, the width of the ground state peak becomes 
narrower with decreasing lattice spacing. We also find that the second peak should be lattice or MEM artifacts since 
its locations varies a lot from lattice spacings and details in this frequency region can not be resolved on the coarsest 
lattices used in this study. 

One always has to compare spectral functions ai T > Tc to those aX T < Tc to study temperature effects. As 
number of correlator data points at higher temperatures is reduced, we study the dependence of output spectral 
functions on the number of data points used in the MEM analysis at T < T^, we use the same number of data points 
below and above Tc to have similar systematic uncertainties and to analyze thermal modifications. Here we restrict 
the default model to have the behavior of the free lattice spectral function. At T = 0.73 Tc we select the data points 
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FIG. 17. Left: the dependence of the output spectral function on the number of correlator data points used in the MEM 
analysis at T = 0.73 Tc- All the points included start at fmin = 4. Af is the step size between the neighboring data points 
selected. For instance, Af = 4 means f = 4, 8, 12, • • • ,48 are used. Right: the fmin (number of data points omitted in the 
short distance) dependence of the output spectral functions at T = 0.73 Tc- "DM" labels the input default model and the other 
lines are the output spectral functions with different values of Af. 



in the temporal direction as to start at fmin = 4 and be separated by a step length of Af . For instance when Af = 2 
wc select data points of f = 4, 6, 8, • • • , 48, in total 13 points. So the number of data points used with Af = 2, 3, 4 
at 0.73 Tc corresponds to the number of data points used at 1.46 2.20 Tc and 2.93 Tc, respectively. We show the 
results for the Va channels in the left plot of Fig. [ITl "DM" labels the input default model and the other lines are 
the output spectral functions with different values of Af . We observe negligible dependences on the number of data 
points used in the interesting frequency region . There are minor changes on the amplitudes of the ground state peak 
but the ground state peak location always remains the same. 

To remove the discretization effects, we normally omit some data points at very small distances. However, it is 
not very certain how many data points should be omitted or up to what value of fmin (the shortest time slice f 
used in the MEM analyses) the physics about the bound states is concerned. Thus we check the dependence of the 
output spectral function on fmin- We vary fmin to be 4, 7, 16, 24 and 36 at 0.73 Tc to check the effects for the same 
default model. The default models are fixed in each channel. The results for the Va channel are shown in the right 
plot of Fig. [ITl "DM" labels the input default model and the other lines are output spectral functions corresponding 
different values of fmin- We observe that the large w (w > 5 GcV) behavior of the output spectral functions, which 
is most sensitive to the small distance part of the correlation function, changes with fmin and in the small w region 
(w ;< 5 GeV) the peak location of the ground state peak stays almost unchanged even with fmin = 36. Thus the fmin 
dependence of the spectral function in the Va channel in the interesting frequency region is very small at T < Tc on 
our finest lattice. 

After exploring the uncertainties that can be learned from correlators at T < Tc, we then move on to study the 
uncertainties of output spectral functions at T > Tc- Following the spirit we have done at T < Tc, we will also 
check the dependences on the number of data points omitted from the short distances. Besides that we will check the 
dependences on the threshold of the continuum through the charm quark mass am. 

We first show the results for the dependence of the output spectral function on the number of correlator data points 
omitted at short distances fmin in Fig. [181 The upper panel shows p{uj)/oj'^ as a function of w at three available 
temperatures above Tc while the lower panel focuses on the transport behavior of the spectral function in the low 
frequency region and has p{uj)/{ujT) as function of lu/T. The default models ("DM") are the same in the whole 
frequency region at each temperature. "DM" is provided by a rescalcd free lattice spectral function and an additional 
transport peak. "DM" is also the same as "DM1" in Fig. [14] at each temperature. Note that the transport part of 
"DM" is parameterized in the Breit-Wigner form with 2ttTD = 3.6 and M = 1.8 GeV. As seen from the upper plot 
of Fig. 1181 at 1.46 Tc, from fmin = 4 to fmin = 7 and 10, the peak location of the ground state peak seems to move a 
little further to larger energy while at both 2.20 Tc and 2.93 Tc the output spectral functions show negligible changes 
due to the variation of fmin = 4, 6, and 7. In the lower panel of Fig. [T8]the very low frequency behavior of the spectral 
function is shown. Note that "DM" in the current frequency region is also fixed at each temperature. At all three 
temperatures the output transport peaks show minor dependences on fmin, which indicates that the information of 
the transport peak is mainly enclosed in the large distance part of the correlation function. 
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FIG. 18. The fmin (number of data points omitted from the short distance) dependence of output spectral functions in the Vu 
channel at T > Tc- The upper plot shows the behavior of p{uj)/uj^ as a function of uj while the lower plot shows the transport 
behavior of p(u])/{u)T) as a function of lo/T which corresponds to the divergent parts in the upper plot at the corresponding 
temperatures. 



Because of the insensitivity of MEM on the very large uj behavior of the spectral function, as we observed from, 
e.g. the left panels of Fig. [HI the outputs always reproduce the very large w behavior of the input default models, 
which in our case normally is the free lattice spectral function multiplied by a certain constant to reproduce the value 
of G{fmin)- However, we do not really know the exact behavior of the large w part as well as the onset point of the 
continuum. We thus check the effects caused by different quark masses in the default models. The different quark 
masses am have an effect on the threshold and the structure of the free spectral function. In Fig. [TH] we show the 
dependence of the output spectral function in the Vu channel on the quark mass am ai T > Tc- The upper panel 
of Fig. [Tni shows the large uj behavior of the spectral function and the lower panel highlights the transport peak 
part. Here we test with free lattice spectral functions having am = 0.06 ("DM1"), am = 0.04 ("DM2") and am = 
("DM3") such that Gj:,m{t ~ 4, r)/G'(f ~ 4,T) ~ 1. Here am = 0.06 is the quark mass obtained from the running 
quark mass on the lattice (see Table |TTI| and "DM1" is the same as "DM1" in Fig. 1181 The rising side of the ground 
state peak starts to be nonzero following the trend of the default model already at 1.46 Tc and the amplitude of the 
ground state peak also changes with different values of am. However, the location of the first peak remains almost the 
same and it is much larger than the ground state peak location at 0.73 Tc- At 2.20 Tc, the output spectral functions 
"spf2" and "spf3" from the default models "DM2" and "DM3" have a small bump structure other than "spfl". At 
2.93 Tc, "spfl" , "spf2" and "spf3" have negligible differences when a; > 3 GeV. As seen from the lower panel of Fig. [TOl 
with decreasing am, in general the transport peak's amplitude becomes smaller and its width becomes larger at all 
the three temperatures. At 1.46 Tc the change of the transport peak of "spfl" is very small, and when going to higher 
temperatures, 2.20 and 2.93 Tc, the deviations become larger, probably as a consequence of the larger differences of 
"spf's in the frequency region of 1 < a; < 7 GeV. 
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FIG. 19. The quark mass am dependence of output spectral functions in tire Va channel at T > Tc. The upper plot shows 
the behavior of p{uj)/uj^ as a function of uj while the lower plot shows the transport behavior of p(u))/{ujT) as function of tj/T 
which corresponds to the divergent parts in the upper plot at the corresponding temperatures. 



In a short summary of this subsection, we have studied the dependence of output spectral functions on the lattice 
spacing and conclude that our finest lattice gives the most reliable results. To better compare the spectral function 
below and above Tc, we used the same number of data points at below and above Tc in the MEM analysis and observed 
negligible difference between the cases with and without the same number of data points used. We checked lattice 
cutoff effects by removing several data points from the short distance and found small dependences of low frequency 
part of spectral function on the short distance of the correlation function. It supports that spectral functions in the 
resonance and transport peak region extracted from MEM using our finest lattice really are physical and show no 
major cutoff dependencies. We also checked the dependences on the threshold of input free lattice spectral function 
and found minor dependences in the resonance part of the spectral function. So one basically sees that the general 
picture in the Va channel is not changed with various different default models and different ways of implementing the 
correlator data. The transport peak in the Va channel has a weak dependence on the threshold of the free lattice 
spectral function at 1.46 and the dependences at T > 1.46Tc become somewhat stronger due to the limited available 
distances. 
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